{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example: Static inverse free-boundary equilibrium calculations\n", "\n", "---\n", "\n", "This example notebook shows how to use FreeGSNKE to solve **static inverse** free-boundary Grad-Shafranov (GS) problems. \n", "\n", "In the **inverse** solve mode we seek to estimate the active poloidal field coil currents using user-defined constraints (e.g. on isoflux, x-point, o-point, and psi values) and plasma current density profiles for a desired equilibrium shape. \n", "\n", "Note that during this solve, currents are **not** found in any specified passive structures. \n", "\n", "Below, we illustrate how to use the solver for both diverted and limited plasma configurations in a **MAST-U-like tokamak** using stored pickle files containing the machine description. These machine description files partially come from the FreeGS repository and are not an exact replica of MAST-U." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The static free-boundary Grad-Shafranov problem\n", "\n", "Here we will outline the **static free-boundary** GS problem that is solved within both the **forward** and **inverse** solvers, though we encourage you to see [Amorisco et al. (2024)](https://pubs.aip.org/aip/pop/article/31/4/042517/3286904/FreeGSNKE-A-Python-based-dynamic-free-boundary) and [Pentland et al. (2024)](https://arxiv.org/abs/2407.12432) for more details. \n", "\n", "Using a cylindrical coordinate system $(R,\\phi,Z)$, the aim is to solve the GS equation:\n", "\n", "$$ \\Delta^* \\psi \\equiv \\left( \\frac{\\partial^2}{\\partial R^2} - \\frac{1}{R} \\frac{\\partial}{\\partial R} + \\frac{\\partial^2}{\\partial Z^2} \\right) \\psi = -\\mu_0 R J_{\\phi}(\\psi, R, Z), \\qquad (R,Z) \\in \\Omega, $$\n", "\n", "for the poloidal flux $\\psi(R,Z)$ (which here has units Weber/$2\\pi$) in the rectangular computational domain $\\Omega$. The flux has contributions from both the plasma and the coils (metals) such that $\\psi = \\psi_p + \\psi_c$. This flux defines the toroidal current density $J_{\\phi} = J_p(\\psi,R,Z) + J_c(R,Z)$, also containing a contribution from both the plasma and coils, respectively. We have the plasma current density (only valid in the core plasma region $\\Omega_p$):\n", "\n", "$$ J_p(\\psi,R,Z) = R \\frac{\\mathrm{d}p}{\\mathrm{d}\\psi} +\\frac{1}{\\mu_0 R} F \\frac{\\mathrm{d} F}{\\mathrm{d} \\psi}, \\qquad (R,Z) \\in \\Omega_p, $$\n", "\n", "where $p(\\psi)$ is the plasma pressure profile and $F(\\psi)$ is the toroidal magnetic field profile. The current density generated by $N$ active coils and passive structures is given by:\n", "\n", "$$ J_c(R,Z) = \\sum_{j=1}^{N} \\frac{I^c_j(R,Z)}{A_j^c}, \\qquad (R,Z) \\in \\Omega, \\quad \\text{where} \\quad I_j^c(R,Z) = \n", "\\begin{cases} \n", " I_j^c & \\text{if } (R,Z) \\in \\Omega_j^c, \\\\ \n", " 0 & \\text{elsewhere}.\n", "\\end{cases}.$$\n", "\n", "This makes use of the current $I^c_j$ in each metal and its cross-sectional area $A^c_j$ (the domain of each metal is denoted by $\\Omega_j^c$).\n", "\n", "To complete the problem, we need the integral (Dirichlet) free-boundary condition:\n", "\n", "$$ \\psi(R,Z) = \\int_{\\Omega} G(R,Z;R',Z') J_{\\phi}(\\psi, R',Z') \\ \\mathrm{d}R' \\mathrm{d}Z', \\qquad (R,Z) \\in \\partial \\Omega, $$\n", "\n", "where $G$ is the (known) Green's function for the elliptic operator above.\n", "\n", "We'll now go through the steps required to solve the **inverse** problem in FreeGSNKE. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create the machine object\n", "\n", "First, we build the machine object from previously created pickle files in the \"machine_configs/MAST-U\" directory. \n", "\n", "FreeGSNKE requires the following paths in order to build the machine:\n", "- `active_coils_path`\n", "- `passive_coils_path`\n", "- `limiter_path`\n", "- `wall_path`\n", "- `magnetic_probe_path` (not required here)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# build machine\n", "from freegsnke import build_machine\n", "tokamak = build_machine.tokamak(\n", " active_coils_path=\"../machine_configs/MAST-U/MAST-U_like_active_coils.pickle\",\n", " passive_coils_path=\"../machine_configs/MAST-U/MAST-U_like_passive_coils.pickle\",\n", " limiter_path=\"../machine_configs/MAST-U/MAST-U_like_limiter.pickle\",\n", " wall_path=\"../machine_configs/MAST-U/MAST-U_like_wall.pickle\",\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the machine\n", "import matplotlib.pyplot as plt\n", "\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n", "plt.tight_layout()\n", "\n", "tokamak.plot(axis=ax1, show=False)\n", "ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "\n", "ax1.grid(alpha=0.5)\n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Instantiate an equilibrium\n", "\n", "We are now ready to build a plasma equilibrium object for our tokamak. This is done using the `freegs4e.Equilibrium` class, which implicitly defines the rectangular domain of the solver as well as the grid resolution.\n", "\n", "`Equilibrium` has sensible defaults, but it is recommended to define the radial and vertical domain of the grid using the `Rmin`, `Rmax`, `Zmin` and `Zmax` parameters, as well as the grid resolution in the radial and vertical directions with the `nx` and `ny` parameters. The grid will be initialised using fourth-order finite differences. Note that the computational grid **must** encompass the domain enclosed by the limiter object (as this is where the plasma will be confined to). Though it does not need to encompass any/all of the active coils or passive structures. \n", "\n", "A tokamak object should be supplied to the `tokamak` parameter to assign the desired machine to the equilibrium.\n", "\n", "If available, an initial guess for the plasma flux $\\psi_p$ (dimensions `nx` x `ny`) can be provided via the `psi` parameter (commented out in the following code). If not, the default initialisation will be used. \n", "\n", "The `eq` object will store a lot of important information and derived quantites once the equilibrium has been calculated (see future notebook on this). " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from freegsnke import equilibrium_update\n", "\n", "eq = equilibrium_update.Equilibrium(\n", " tokamak=tokamak, # provide tokamak object\n", " Rmin=0.1, Rmax=2.0, # radial range\n", " Zmin=-2.2, Zmax=2.2, # vertical range\n", " nx=65, # number of grid points in the radial direction (needs to be of the form (2**n + 1) with n being an integer)\n", " ny=129, # number of grid points in the vertical direction (needs to be of the form (2**n + 1) with n being an integer)\n", " # psi=plasma_psi\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Instantiate a profile object\n", "\n", "We can now instatiate a profile object that contains the chosen parameterisation of the toroidal plasma current density $J_p$ (i.e. on right hand side of the GS equation). We can then set the paramters for the chosen current density profiles. \n", "\n", "The following table indicates which parameterisations are currently available in FreeGSNKE and for which types of equilibrium simulation they currently work for:\n", "\n", "| Name | Static | Evolutive (linear) | Evolutive (nonlinear) |\n", "|-------|--------|------------------|---------------------|\n", "| `ConstrainPaxisIp` | ✅ | ✅ | ✅ |\n", "| `ConstrainBetapIp` | ✅ | ✅ | ✅ |\n", "| `FiestaTopeol` | ✅ | ✅ | ✅ |\n", "| `Lao85` | ✅ | ✅ | ✅ |\n", "| `TensionSpline` | ✅ | ❌ | ❌ |\n", "| `GeneralPprimeFFprime` | ✅ | ❌ | ❌ |\n", "\n", "In this notebook, we will make use of the `ConstrainPaxisIp` (and `ConstrainBetapIp`) profiles (see [Jeon (2015)](https://link.springer.com/article/10.3938/jkps.67.843)). Others will be utilised in later notebooks. If there is a profile parameterisation you require that does not exist, please do create an issue.\n", "\n", "Both `ConstrainPaxisIp` and `ConstrainBetapIp` are parameterised as follows:\n", " $$J_{p}(\\psi, R, Z) = \\lambda\\big[ \\beta_{0} \\frac{R}{R_{0}} \\left( 1-\\tilde{\\psi}^{\\alpha_m} \\right)^{\\alpha_n} + (1-\\beta_{0}) \\frac{R_0}{R} \\left( 1-\\tilde{\\psi}^{\\alpha_m} \\right)^{\\alpha_n} \\big] \\quad (R,Z) \\in \\Omega_p,$$\n", "\n", "where the first term is the pressure profile and the second is the toroidal current profile. Here, $\\tilde{\\psi}$ denotes the normalised flux:\n", "$$ \\tilde{\\psi} = \\frac{\\psi - \\psi_a}{\\psi_b - \\psi_a}, $$\n", "where $\\psi_a$ and $\\psi_b$ are the values of the flux on the magnetic axis and plasma boundary, respectively. \n", "\n", "The parameters required to define this particular profile are:\n", "- `Ip` (total plasma current).\n", "- `fvac` ($R B_{tor}$, vacuum toroidal field strength).\n", "- `alpha_m`>0, and `alpha_n`>0 (that define the shape/peakedness of the profiles).\n", "- If `ConstrainPaxisIp` is used, then `paxis` (pressure on the magnetic axis) is required.\n", "- If `ConstrainBetapIp` is used, then `betap` (proxy of the poloidal beta) is required.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The values of $\\lambda$ and $\\beta_0$ are found using the above parameters as constraints ($R_0$ is a fixed scaling constant) in the following:\n", "\n", "- For `ConstrainPaxisIp`, we can re-arrange the following equations to solve for the unknowns:\n", "\n", "$$ p_{\\text{axis}} = \\lambda \\beta_{0} \\frac{R}{R_{0}} \\int^{\\psi_b}_{\\psi_a} \\left( 1-\\tilde{\\psi}^{\\alpha_m} \\right)^{\\alpha_n} \\mathrm{d} \\tilde{\\psi} $$\n", "\n", "and\n", "\n", "$$ I_p = \\int^{Z_{\\text{max}}}_{Z_{\\text{min}}} \\int^{R_{\\text{max}}}_{R_{\\text{min}}} J_p(\\psi, R, Z) \\ \\mathrm{d}R \\mathrm{d}Z. $$\n", "\n", "\n", "- For `ConstrainBetapIp`, we can instead re-arrange and solve the following:\n", "\n", "$$ \\beta_p = \\frac{8 \\pi}{\\mu_0 I_p^2} \\int^{Z_{\\text{max}}}_{Z_{\\text{min}}} \\int^{R_{\\text{max}}}_{R_{\\text{min}}} p(\\psi) \\ \\mathrm{d}R \\mathrm{d}Z. $$\n", "\n", "and\n", "\n", "$$I_p = \\int^{Z_{\\text{max}}}_{Z_{\\text{min}}} \\int^{R_{\\text{max}}}_{R_{\\text{min}}} J_p(\\psi, R, Z) \\ \\mathrm{d}R \\mathrm{d}Z. $$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In what follows, we use `ConstrainPaxisIp`. Note that the equilibrium (`eq`) object is passed to the profile to inform calculations relating to the machine description." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# initialise the profiles\n", "from freegsnke.jtor_update import ConstrainPaxisIp\n", "profiles = ConstrainPaxisIp(\n", " eq=eq, # equilibrium object\n", " paxis=8e3, # pressure on the magnetic axis\n", " Ip=6e5, # plasma current\n", " fvac=0.5, # fvac = R B_{tor}\n", " alpha_m=1.8, # profile function parameter\n", " alpha_n=1.2 # profile function parameter\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load the static nonlinear solver\n", "\n", "We can now load FreeGSNKE's static solver object. The equilibrium is used to inform the solver of the computational domain and of the tokamak properties. The solver below can be used for **both** inverse and forward solves. \n", "\n", "Note: It's not necessary to instantiate a new solver when aiming to use it on new or different equilibria, as long as the integration domain, mesh grid, and tokamak are consistent across solves. \n", "\n", "Note: A random seed can also be set here. This comes into play for tricky forward/inverse solves where the internal Newton-Krylov solver gets stuck and needs to search in random Krylov direction. See code documentation for further details." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from freegsnke import GSstaticsolver\n", "GSStaticSolver = GSstaticsolver.NKGSsolver(\n", " eq=eq, # eq object\n", " # seed=42, # random seed\n", " ) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fixed coil currents\n", "\n", "Firstly, we set the values for any currents in specific active poloidal field coils that we may know. The inverse solver will not allow the current in these coils to vary if the `control` parameter is set to `False` (i.e. it will be excluded during the optimisation). \n", "\n", "Note: any passive structures in the tokamak automatically have their control parameter set to False and are therefore not included in an inverse solve. \n", "\n", "As an example, we will fix the `Solenoid` current and seek a solution in which this value is fixed, rather than estimated by the inverse solve." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "eq.tokamak.set_coil_current('Solenoid', 5000)\n", "eq.tokamak['Solenoid'].control = False # ensures the current in the Solenoid is fixed" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Constraints\n", "\n", "Recall that in the **inverse** solve mode we seek to **estimate the active poloidal field coil currents** using user-defined **constraints** (e.g. on isoflux, null points, and $\\psi$ values) together with prescribed plasma current density profiles, in order to obtain a desired equilibrium shape.\n", "\n", "FreeGSNKE uses a `constrain` object, which accepts the following types of constraints:\n", "\n", "- **Null points (`null_points`)** \n", " A null point $(R_X, Z_X)$ (X-point or O-point) is defined by vanishing poloidal field, i.e. we can impose\n", " $$\n", " \\nabla \\psi(R_X, Z_X) = \\mathbf{0}.\n", " $$\n", " In particular, an X-point is a saddle point of $\\psi$ while an O-point (magnetic axis) is a local extremum (typically a minimum) of $\\psi$. \n", " The above is equivalent to both:\n", " $$\n", " \\frac{\\partial \\psi}{\\partial R}(R_X, Z_X) = \\frac{\\partial \\psi}{\\partial Z}(R_X, Z_X) = 0,\n", " $$\n", " or\n", " $$\n", " B_r (R_X, Z_X) = B_z(R_X, Z_X) = 0.\n", " $$\n", " This means that for each null point considered, two constraints are imposed internall within the inverse solver. \n", " \n", "\n", "- **Isoflux sets (`isoflux_set`)** \n", " An isoflux set is a set of points $\\{(R_i, Z_i)\\}_{i=1}^k$ that lie on the same flux surface, i.e. we can impose\n", " $$\n", " \\psi(R_1, Z_1) = \\psi(R_2, Z_2) = \\cdots = \\psi(R_k, Z_k).\n", " $$\n", " Internally, the inverse solver enforces each constraint such that\n", " $$\n", " \\psi(R_i, Z_i) - \\psi(R_j, Z_j) = 0,\n", " \\qquad i \\neq j,\n", " $$\n", " meaning that for $k$ points in an isoflux set, there are $k(k-1)/2$ constraints. \n", "\n", "- **Coil current limits (`coil_current_limits`)** \n", " If we wnt to enforce bounds on the coil currents (for machine safety purposes), we can impose\n", " $$\n", " I_m^{\\min} \\le I_m \\le I_m^{\\max},\n", " \\qquad m = 1, \\dots, N.\n", " $$\n", "\n", "**Note**: see the next notebooks for some more advanced constraints!\n", "\n", "At least one constraint (preferably more) is required to carry out an inverse solve. Although it may be tempting to utilise all of the above constraints, one must be careful not to **overconstrain** the inverse problem, else constraints will fight one another and the solver will not converge to an equilibrium. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the following, we start simple and specify two X-point locations (to obtain a double-null plasma), a few isoflux locations, and some coil current limits. The isoflux constraints define the core plasma boundary and part of the divertor legs." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np \n", "from freegsnke.inverse import Inverse_optimizer\n", "\n", "Rx = 0.6 # X-point radius\n", "Zx = 1.1 # X-point height\n", "Rout = 1.4 # outboard midplane radius\n", "Rin = 0.34 # inboard midplane radius\n", "\n", "# set desired null_points locations (this can include X-point and O-point locations)\n", "null_points = [[Rx, Rx], [Zx, -Zx]]\n", "\n", "# set desired isoflux constraints with format \n", "# isoflux_set = [isoflux_0, isoflux_1 ... ] \n", "# with each isoflux_i = [R_coords, Z_coords]\n", "isoflux_set = np.array([\n", " [[Rx, Rx, Rin, Rout, 1.0, 1.0, .8,.8], [Zx, -Zx, 0.,0., 2.0, -2.0, 1.62, -1.62]]\n", " ])\n", "\n", "# set the coil current limits (upper and lower)\n", "# coil ordering in this case: PX, D1, D2, D3, Dp, D5, D6, D7, P4, P5, P6\n", "coil_current_limits = [\n", " [5e3, 9e3, 9e3, 7e3, 7e3, 5e3, 4e3, 5e3, 0.0, 0.0, None],\n", " [-5e3, -9e3, -9e3, -7e3, -7e3, -5e3, -4e3, -5e3, -10e3, -10e3, None]\n", "]\n", " \n", "# instantiate the freegsnke constrain object\n", "constrain = Inverse_optimizer(\n", " null_points=null_points,\n", " isoflux_set=isoflux_set,\n", " coil_current_limits=coil_current_limits\n", ")\n", "\n", "# if you find coil limits are being violated or an undesireable solution is being produced,\n", "# you can try increasing the penalty factor for violating the coil limits\n", "# (here we just set it to its default value of 1e5)\n", "constrain.mu_coils = 1e5" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The linear system of constraints\n", "\n", "During an inverse solve, a minmisation problem involving the responses, changes in coil currents, and constraints, is repeatedly solved:\n", "\n", "$$ \\min_{x} \\| A\\vec{x} - \\vec{b}\\|^2 + \\| \\gamma \\vec{x} \\|^2, $$\n", "\n", "where\n", "- $A$ is the fixed response matrix (that determines how a change in coil currents $x$, affects constraint values). \n", "- $\\vec{x} = \\Delta \\vec{I}^c$ is the step change in active coil currents required to match the constraints. \n", "- $\\vec{b}$ is the vector of constraint values being enforced. \n", "- $\\gamma > 0$ is the regularisation parameter/vector.\n", "\n", "We solve for $\\vec{x}$ using a gradient-based optimiser. \n", "\n", "[Song et al. (2024)](https://www.mdpi.com/2571-6182/7/4/45) provide a nice overview of the inverse problem.\n", "\n", "The coil limits are handled in a slightly different manner because they are box constraints that cannot be encoded in the minimisation problem. Instead, FreeGSNKE adds additional constraints to the solver that bound the solution $\\Delta \\vec{I}^c$. These bounds have a tolerance (slack) which allows the coils to violate their limits, but penalises violations; this allows the coil limits to be violated 'on the path' to a solution that respects the prescribed limits.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The inverse solve\n", "\n", "The following cell will execute the solve. Since a `constrain` object is provided, this is interpreted as a call to the inverse solver, if `constrain=None`, then the forward solver will be called (see next notebook). \n", "\n", "The `target_relative_tolerance` is the maximum relative error on the plasma flux function allowed for convergence and `target_relative_psit_update` ensures that the relative update to the plasma flux (caused by the update in the control currents) is lower than this target value. Both are required to be met for the inverse problem to be considered successfully solved.\n", "\n", "The `verbose=True` option will provide an indication of the progression of the solve. \n", "\n", "The `l2_reg` parameter defines the Tikonov regularisation (i.e. $\\gamma$ in the cell above) used by the optmiser. This can be set as a scalar or as a vector (equal to the number of coil currents being solved for). Larger values force coil currents to stay closer to their original values while lower values will encourage more freedom. This can be useful in particular for vertical control coils (e.g. the P6 coil in MAST-U), in which we don't want the coil current to \"jump around\" during optimisation and cause vertical instability issues. \n", "\n", "\n", "The solver steps are (roughly):\n", "1. Solve the linear system to find initial coil currents that approximately satisfy the constraints for the initial equilibrium. \n", "2. While tolerance is not met:\n", " - use the coil currents to solve the forward GS problem (with Newton-Krylov iterations).\n", " - solve the linear system to update coil currents to satisfy constraints for current equilibrium (from forward solve). \n", " - check tolerances. \n", "\n", "Given that there may be more or less constraints than unknown parameters (i.e. coil currents), the inverse problem may be over- or under-constrained. This means that there may be zero, one, or many solutions to the problem. It is often a bit of an art when solving inverse problems so the best strategy is to start with some basic isoflux constraints before trying to add null points and others." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# solve!\n", "GSStaticSolver.solve(eq=eq, \n", " profiles=profiles, \n", " constrain=constrain, \n", " target_relative_tolerance=1e-6,\n", " target_relative_psit_update=1e-3,\n", " verbose=True, # print output\n", " l2_reg=np.array([1e-12]*10+[1e-6]), \n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following plots show how to display:\n", "1. the tokamak with:\n", " - active coil filaments (rectangles with blue interior)\n", " - passive structures (blue circles if defined as filaments or thin black outline/grey interior if defined as parallelograms)\n", " - limiter/wall (solid black)\n", "2. the tokamak + the equilibrium with:\n", " - separatrix/last closed flux surface (solid red lines)\n", " - poloidal flux (yellow/green/blue contours, colours indicates magnitude)\n", " - X-points (red circles)\n", " - O-points (green crosses)\n", "3. the tokamak + the equilibrum + the contraints with:\n", " - null-point constraints (purple Y-shaped crosses)\n", " - isoflux contour constraints (usual crosses)\n", "\n", "Note: Setting 'show=True' can toggle the legend on/off." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 8), dpi=80)\n", "\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.set_aspect('equal')\n", "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "\n", "ax2.grid(zorder=0, alpha=0.75)\n", "ax2.set_aspect('equal')\n", "eq.tokamak.plot(axis=ax2,show=False) # plots the active coils and passive structures\n", "ax2.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", "eq.plot(axis=ax2,show=False) # plots the equilibrium\n", "ax2.set_xlim(0.1, 2.15)\n", "ax2.set_ylim(-2.25, 2.25)\n", "\n", "\n", "ax3.grid(zorder=0, alpha=0.75)\n", "ax3.set_aspect('equal')\n", "eq.tokamak.plot(axis=ax3,show=False) # plots the active coils and passive structures\n", "ax3.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", "eq.plot(axis=ax3,show=False) # plots the equilibrium\n", "constrain.plot(axis=ax3, show=True) # plots the contraints\n", "ax3.set_xlim(0.1, 2.15)\n", "ax3.set_ylim(-2.25, 2.25)\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A solve call will modify the equilibrium object in place. That means that certain quantities within the object will be updated as a result of the solve. \n", "\n", "Various different quantities and functions can be accessed via the 'eq' and 'profile' objects. For example:\n", "- the total flux can be accessed with `eq.psi()`.\n", "- the plasma flux with `eq.plasma_psi`.\n", "- the active coil + passive structure flux with `eq.tokamak_psi`.\n", "- (Total flux = plasma flux + coil flux)\n", "\n", "Explore `eq.` to see more (also `profiles.`, e.g. the plasma current distribution over the domain can be found with `profiles.jtor`). Also see the `example3` notebook for how to access many different equilibrium-derived quantities of interest. \n", "\n", "The set of optimised coil currents can be retrieved using `eq.tokamak.getCurrents()` having been assigned to the equilibrium object during the inverse solve.\n", "\n", "The following lines will save the calculated currents to a pickle file (we will use these in future notebooks). " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "inverse_currents = eq.tokamak.getCurrents()\n", "inverse_currents_arr = eq.tokamak.getCurrentsVec()[0:12]\n", "\n", "# save coil currents to file\n", "import pickle\n", "with open('data/simple_diverted_currents_PaxisIp.pk', 'wb') as f:\n", " pickle.dump(obj=inverse_currents, file=f)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we can check that the currents we have found respect the coil limits that we prescribed earlier:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# storage\n", "coil_names = []\n", "active_currents = []\n", "min_limits = []\n", "max_limits = []\n", "\n", "for coil_name, coil_current, ul, ll in zip(\n", " eq.tokamak.coils_dict,\n", " eq.tokamak.getCurrentsVec()[:12],\n", " [None] + coil_current_limits[0], # upper limits\n", " [None] + coil_current_limits[1], # lower limits\n", "):\n", " # skip solenoid (no limits)\n", " if ul is None or ll is None:\n", " continue\n", "\n", " coil_names.append(coil_name)\n", " active_currents.append(coil_current)\n", " min_limits.append(ll)\n", " max_limits.append(ul)\n", " \n", " if (coil_current >= ll) & (coil_current <= ul):\n", " within_limit = True\n", " else:\n", " within_limit = False\n", "\n", " # print\n", " print(f\"{coil_name}: {coil_current:.2f} --> [{ll},{ul}] [A] --> limits met = {within_limit}\")\n", "\n", "# convert to numpy\n", "active_currents = np.array(active_currents)\n", "min_limits = np.array(min_limits)\n", "max_limits = np.array(max_limits)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also visualise the currents with respect to their limits (normalising between -1 and 1). Note the Solenoid and P6 are missing as they have no limits in this case. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# VISUALISE HOW CLOSE THEY ARE TO THE LIMITS (normalised)\n", "\n", "# normalized to [-1, 1] scale (relative to min/max)\n", "norm_current = 2 * (active_currents - min_limits) / (max_limits - min_limits) - 1\n", "\n", "# plot\n", "plt.figure(figsize=(12, 5), dpi=70)\n", "x = np.arange(len(coil_names))\n", "\n", "# line plot\n", "plt.plot(x, norm_current, marker='x', linestyle='-', color='red', label=\"Current\")\n", "plt.axhline(-1, color='k', linestyle='-', linewidth=1, label=\"Min. limit\")\n", "plt.axhline(1, color='k', linestyle='-', linewidth=1, label=\"Max. limit\")\n", "\n", "# formatting\n", "plt.xticks(x, coil_names, rotation=45, ha='right')\n", "plt.ylabel('Coil name [Amps]')\n", "plt.ylabel('Normalised current')\n", "plt.title('Coil current proximity to limits (normalised)')\n", "plt.grid(True, linestyle='--', alpha=0.6)\n", "plt.legend()\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Second inverse solve: limited plasma example\n", "\n", "Below we carry out an inverse solve seeking coil current values for a limited plasma configuration (rather than a diverted one). \n", "\n", "In a limiter configuration the plasma \"touches\" the limiter of the tokamak and is confined by the solid structures of the vessel. The last closed flux surface (LCFS) is the closed contour that is furthest from the magnetic axis that just barely touches (and is tangent to) the limiter." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# a new eq object resets the coil currents and equilibrium\n", "eq = equilibrium_update.Equilibrium(\n", " tokamak=tokamak, # provide tokamak object\n", " Rmin=0.1, Rmax=2.0, # radial range\n", " Zmin=-2.2, Zmax=2.2, # vertical range\n", " nx=65, # number of grid points in the radial direction (needs to be of the form (2**n + 1) with n being an integer)\n", " ny=129, # number of grid points in the vertical direction (needs to be of the form (2**n + 1) with n being an integer)\n", ")\n", "\n", "# initialise the profiles\n", "profiles = ConstrainPaxisIp(\n", " eq=eq, # equilibrium object\n", " paxis=6e3, # pressure on the magnetic axis <-- changed\n", " Ip=4e5, # plasma current <-- changed\n", " fvac=0.5, # fvac = R B_{tor}\n", " alpha_m=1.8, # profile function parameter\n", " alpha_n=1.2 # profile function parameter\n", ")\n", "\n", "# first we specify some alternative constraints\n", "Rout = 1.4 # outboard midplane radius\n", "Rin = 0.255 # inboard midplane radius (less than the wall radius to make it limited)\n", "\n", "# locations of X-points\n", "Rx = 0.45\n", "Zx = 1.18\n", "null_points = [[Rx, Rx], [Zx, -Zx]]\n", "\n", "# isoflux constraints\n", "isoflux_set = [\n", " [Rx, Rx, Rout, Rin, Rin, Rin, Rin, Rin, .75, .75, .85, .85],\n", " [Zx, -Zx, 0, 0, .1, -.1, .2, -.2, 1.6, -1.6, 1.7, -1.7]\n", " ]\n", "\n", "# let's assume we're seeking an equilibrium with no solenoid current\n", "eq.tokamak.set_coil_current('Solenoid', 0)\n", "eq.tokamak['Solenoid'].control = False # fixes the current\n", "\n", "# also make it perfectly up/down symmetric (as P6 is wired in anti-series vs. the other coils)\n", "eq.tokamak.set_coil_current('P6', 0)\n", "eq.tokamak['P6'].control = False # fixes the current\n", "\n", "# pass the magnetic constraints to a new constrain object\n", "constrain = Inverse_optimizer(\n", " null_points=null_points,\n", " isoflux_set=isoflux_set\n", " )\n", "\n", "# carry out the solve\n", "GSStaticSolver.solve(eq=eq, \n", " profiles=profiles, \n", " constrain=constrain, \n", " target_relative_tolerance=1e-6,\n", " target_relative_psit_update=1e-3,\n", " verbose=True, # print output\n", " l2_reg=np.array([1e-9]*10),\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# save the currents for later use\n", "inverse_current_values = eq.tokamak.getCurrents()\n", "\n", "# save coil currents to file\n", "with open('data/simple_limited_currents_PaxisIp.pk', 'wb') as f:\n", " pickle.dump(obj=inverse_current_values, file=f)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the resulting equilbria \n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n", "ax1.grid(True, which='both')\n", "eq.plot(axis=ax1, show=False)\n", "eq.tokamak.plot(axis=ax1, show=False)\n", "constrain.plot(axis=ax1,show=True)\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": ".venv", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.16" } }, "nbformat": 4, "nbformat_minor": 4 }